The smuthi.linearsystem.particlecoupling package

linearsystem.particlecoupling

Routines for multiple scattering. Some modules contain functions to explicitly compute the coupling matrix entries. Others contains functions for the preparation of lookup tables that are used to approximate the coupling matrices by interoplation. The Cython folder has functions to calculate the direct coupling block between two particles in the same layer using fast C routines.

particlecoupling.direct_coupling

This module contains functions to compute the direct (i.e., not layer mediated) particle coupling coefficients.

smuthi.linearsystem.particlecoupling.direct_coupling.direct_coupling_block(vacuum_wavelength, receiving_particle, emitting_particle, layer_system)
Direct particle coupling matrix \(W\) for two particles that do not have intersecting circumscribing spheres.

This routine is explicit.

To reduce computation time, this routine relies on two internal accelerations. First, in most cases the number of unique maximum multipole indicies, \((\tau, l_{max}, m_{max})\), is much less than the number of unique particles. Therefore, all calculations that depend only on multipole indicies are stored in an intermediate hash table. Second, Cython acceleration is used by default to leverage fast looping. If the Cython files are not supported, this routine will fall back on equivalent Python looping.

Cython acceleration can be between 10-1,000x faster compared to the Python equivalent. Speed variability depends on the number of unique multipoles indicies, the size of the largest multipole order, and if particles share the same z coordinate.

Parameters:
Returns:

Direct coupling matrix block as numpy array.
smuthi.linearsystem.particlecoupling.direct_coupling.direct_coupling_block_2D_from_hash_table(blocksize1, blocksize2, w, sph, k, dx, dy, dz, lmax1, lmax2, mmax1, mmax2, a5leg_hash_table, b5leg_hash_table, threaded=False)
Subroutine to calculate the direct coupling between two particles
that are in the same layer and have the same z coordinate. This subroutine is called internally by direct_coupling_block. If Cython is enabled then this subroutine is Cython accelerated. Otherwise, the Python equivalent subroutine is used.
Parameters:
  • blocksize1 (int) – Number of columns in the direct coupling block
  • blocksize2 (int) – Number of rows in the direct coupling block
  • w (ndarray) – Zero initialized direct coupling matrix block as numpy array
  • sph (ndarray) – Zero initialized Hankel function matrix as numpy array
  • k (complex) – Wavenumber in the shared media
  • dx (float) – x-coordinate of the distance between the two particles
  • dy (float) – y-coordinate of the distance between the two particles
  • dz (float) – z-coordinate of the distance between the two particles
  • lmax1 (int) – Largest polar quantum number of the recieving particle
  • lmax2 (int) – Largest polar quantum number of the emitting particle
  • mmax1 (int) – Largest azimuthal quantum number of the recieving particle
  • mmax2 (int) – Largest azimuthal quantum number of the emitting particle
  • a5leg_array (ndarray) – Hash table of the elements in \(A\) not dependent on \(/phi\) or \(kd\)
  • b5leg_array (ndarray) – Hash table of the elements in \(B\) not dependent on \(/phi\) or \(kd\)
  • threaded (bool) – Flag to enable multithreading (valid only for Cython where the Gil can be released). Currently hard coded to False.
Returns:

Direct coupling matrix block as numpy array.
Return type:

w (ndarray)
smuthi.linearsystem.particlecoupling.direct_coupling.direct_coupling_block_3D_from_hash_table(blocksize1, blocksize2, w, sph, k, dx, dy, dz, lmax1, lmax2, mmax1, mmax2, a5_hash_table, b5_hash_table, threaded=False)
Subroutine to calculate the direct coupling between two particles
that are in the same layer and do not have the same z coordinate. This subroutine is called internally by direct_coupling_block. If Cython is enabled then this subroutine is Cython accelerated. Otherwise, the Python equivalent subroutine is used.
Parameters:
  • blocksize1 (int) – Number of columns in the direct coupling block
  • blocksize2 (int) – Number of rows in the direct coupling block
  • w (ndarray) – Zero initialized direct coupling matrix block as numpy array
  • sph (ndarray) – Zero initialized Hankel function matrix as numpy array
  • k (complex) – Wavenumber in the shared media
  • dx (float) – x-coordinate of the distance between the two particles
  • dy (float) – y-coordinate of the distance between the two particles
  • dz (float) – z-coordinate of the distance between the two particles
  • lmax1 (int) – Largest polar quantum number of the recieving particle
  • lmax2 (int) – Largest polar quantum number of the emitting particle
  • mmax1 (int) – Largest azimuthal quantum number of the recieving particle
  • mmax2 (int) – Largest azimuthal quantum number of the emitting particle
  • a5_array (ndarray) – Hash table of the elements in \(A\) not dependent on \(/theta\), \(/phi\) or \(kd\)
  • b5_array (ndarray) – Hash table of the elements in \(B\) not dependent on \(/theta\), \(/phi\) or \(kd\)
  • threaded (bool) – Flag to enable multithreading (valid only for Cython where the Gil can be released). Currently hard coded to False.
Returns:

Direct coupling matrix block as numpy array.
Return type:

w (ndarray)
smuthi.linearsystem.particlecoupling.direct_coupling.direct_coupling_block_pvwf_mediated(vacuum_wavelength, receiving_particle, emitting_particle, layer_system, k_parallel, alpha=None, beta=None)

Direct particle coupling matrix \(W\) for two particles (via plane vector wave functions). For details, see: Dominik Theobald et al., Phys. Rev. A 96, 033822, DOI: 10.1103/PhysRevA.96.033822 or arXiv:1708.04808

The plane wave coupling is performed in a rotated coordinate system, which must be chosen such that both particles can be separated by a plane that is parallel to the xy-plane (such that the emitting particle is entirely above that plane and the receiving particle is entirely below that plane).

Two angles (alpha and beta) are required to specify the active rotation into that coordinate system, i.e., the rotation which rotates the particle locations such that the abovementioned condition is fulfilled.

If the angle arguments are omitted, Smuthi tries to estimate a suitable separating plane.

Parameters:
  • vacuum_wavelength (float) – Vacuum wavelength \(\lambda\) (length unit)
  • receiving_particle (smuthi.particles.Particle) – Particle that receives the scattered field
  • emitting_particle (smuthi.particles.Particle) – Particle that emits the scattered field
  • layer_system (smuthi.layers.LayerSystem) – Stratified medium in which the coupling takes place
  • k_parallel (numpy.array) – In-plane wavenumber for plane wave expansion
  • alpha (float) – First Euler angle, rotation around z-axis, in rad
  • beta (float) – Second Euler angle, rotation around y’-axis in rad
Returns:

Direct coupling matrix block (numpy array).
smuthi.linearsystem.particlecoupling.direct_coupling.direct_coupling_matrix(vacuum_wavelength, particle_list, layer_system)

Return the direct particle coupling matrix W for a particle collection in a layered medium.

Parameters:
  • vacuum_wavelength (float) – Wavelength in length unit
  • (list of smuthi.particles.Particle obejcts (particle_list) – Scattering particles
  • layer_system (smuthi.layers.LayerSystem) – The stratified medium
Returns:

Ensemble coupling matrix as numpy array.
smuthi.linearsystem.particlecoupling.direct_coupling.get_separating_plane(receiving_particle, emitting_particle)

Estimate the orientation of a plane that separates two particles. Such a plane always exists for convex particles.

For some special cases (pairs of spheroids, pairs of non-rotated cylinders) a separating plane is constructed.

For all other cases, the orientation is chosen along the center-to-center vector. Note that this choice is not guaranteed to yield a correct PVWF coupling result.

Args: receiving_particle (smuthi.particles.Particle): Receiving particle emitting_particle (smuthi.particles.Particle): Emitting particle
Retruns:

Tuple containing Euler angles for the active rotation into a frame such that the emitting particle is above some z-plane and the receiving particle is below that plane:

  • first rotation Euler angle alpha (float)
  • second rotation Euler angle beta (float)
smuthi.linearsystem.particlecoupling.direct_coupling.spheroids_closest_points(ab_halfaxis1, c_halfaxis1, center1, orientation1, ab_halfaxis2, c_halfaxis2, center2, orientation2)

Computation of the two closest points of two adjacent spheroids. For details, see: Stephen B. Pope, Algorithms for Ellipsoids, Sibley School of Mechanical & Aerospace Engineering, Cornell University, Ithaca, New York, February 2008

Parameters:
  • ab_halfaxis1 (float) – Half axis orthogonal to symmetry axis of spheroid 1
  • c_halfaxis1 (float) – Half axis parallel to symmetry axis of spheroid 1
  • center1 (numpy.array) – Center coordinates of spheroid 1
  • orientation1 (numpy.array) – Orientation angles of spheroid 1
  • ab_halfaxis2 (float) – Half axis orthogonal to symmetry axis of spheroid 2
  • c_halfaxis2 (float) – Half axis parallel to symmetry axis of spheroid 2
  • center2 (numpy.array) – Center coordinates of spheroid 2
  • orientation2 (numpy.array) – Orientation angles of spheroid 2
Retruns:
Tuple containing:
  • closest point on first particle (numpy.array)
  • closest point on second particle (numpy.array)
  • first rotation Euler angle alpha (float)
  • second rotation Euler angle beta (float)

particlecoupling.layer_mediated_coupling

This module contains functions to compute the layer mediated particle coupling coefficients.

smuthi.linearsystem.particlecoupling.layer_mediated_coupling.eval_BeLBe
smuthi.linearsystem.particlecoupling.layer_mediated_coupling.g_function(vacuum_wavelength, receiving_particle, emitting_particle, layer_system, k_parallel)
This function returns the function g_n(k_
ho) as defined in equation (16) of the paper “A quick way to approximate

a Sommerfeld-Weyl_type Sommerfeld integral” by Chew (1988) for the Sommerfeld integral for the layer-mediated particle coupling (compare Amos Egel’s dissertation, equation (3.46)). The purpose of this function is to allow for stationary phase approximation in the case of large lateral distances.

The use of this function is constrained to the important special case of particles above a (possibly layered) substrate.

Args:

vacuum_wavelength (float): Vacuum wavelength \(\lambda\) (length unit) receiving_particle (smuthi.particles.Particle): Particle that receives the scattered field (must be in top

layer)
emitting_particle (smuthi.particles.Particle): Particle that emits the scattered field (must be in top
layer)

layer_system (smuthi.layers.LayerSystem): Stratified medium in which the coupling takes place k_parallel (float): In-plane wavenumber

Returns:
a tuple with the following items: - g: a matrix of g_function as entries (dimension is blocksize x blocksize) - rho: radial distance between particles - delta_z: z-argument of e^(ik_z z) term - bessel_order: matrix of positive numbers to be used as the order of Hankel functions
smuthi.linearsystem.particlecoupling.layer_mediated_coupling.layer_mediated_coupling_block(vacuum_wavelength, receiving_particle, emitting_particle, layer_system, k_parallel='default', show_integrand=False)

Layer-system mediated particle coupling matrix \(W^R\) for two particles. This routine is explicit, but slow.

Parameters:
  • vacuum_wavelength (float) – Vacuum wavelength \(\lambda\) (length unit)
  • receiving_particle (smuthi.particles.Particle) – Particle that receives the scattered field
  • emitting_particle (smuthi.particles.Particle) – Particle that emits the scattered field
  • layer_system (smuthi.layers.LayerSystem) – Stratified medium in which the coupling takes place
  • k_parallel (numpy ndarray) – In-plane wavenumbers for Sommerfeld integral If ‘default’, use smuthi.fields.default_Sommerfeld_k_parallel_array
  • show_integrand (bool) – If True, the norm of the integrand is plotted.
Returns:

Layer mediated coupling matrix block as numpy array.
smuthi.linearsystem.particlecoupling.layer_mediated_coupling.layer_mediated_coupling_block_stat_phase_approx(vacuum_wavelength, receiving_particle, emitting_particle, layer_system)

Compute the layer mediated coupling coefficients by means of the stationary phase approximation, as presented in “A quick way to approximate a Sommerfeld-Weyl_type Sommerfeld integral” by W.C. Chew (1988).

The stationary phase approximation is expected to yield good results for particles with a large lateral distance.

Note: This function assumes that both particles (emitter and receiver) are located in the top layer of the layered medium. For other cases, this function does not apply. ****************************************************************************************************************

Parameters:
Returns:

Stationary phase approximation for the layer mediated coupling matrix block as numpy array.
smuthi.linearsystem.particlecoupling.layer_mediated_coupling.layer_mediated_coupling_matrix(vacuum_wavelength, particle_list, layer_system, k_parallel='default')

Layer system mediated particle coupling matrix W^R for a particle collection in a layered medium.

Parameters:
  • vacuum_wavelength (float) – Wavelength in length unit
  • (list of smuthi.particles.Particle obejcts (particle_list) – Scattering particles
  • layer_system (smuthi.layers.LayerSystem) – The stratified medium
  • k_parallel (numpy.ndarray or str) – In-plane wavenumber for Sommerfeld integrals. If ‘default’, smuthi.fields.default_Sommerfeld_k_parallel_array
Returns:

Ensemble coupling matrix as numpy array.
smuthi.linearsystem.particlecoupling.layer_mediated_coupling.numba_trapz

particlecoupling.prepare_lookup

This module contains functionality to prepare lookups for the particle coupling coefficients, which allows to efficiently treat large numbers of particles.

smuthi.linearsystem.particlecoupling.prepare_lookup.radial_coupling_lookup_table(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None)

Prepare Sommerfeld integral lookup table to allow for a fast calculation of the coupling matrix by interpolation. This function is called when all particles are on the same z-position.

Parameters:
  • vacuum_wavelength (float) – Vacuum wavelength in length units
  • particle_list (list) – List of particle objects
  • layer_system (smuthi.layers.LayerSystem) – Stratified medium
  • k_parallel (numpy.ndarray or str) – In-plane wavenumber for Sommerfeld integrals. If ‘default’, smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float) – Spatial resolution of lookup table in length units. (default: vacuum_wavelength / 100) Smaller means more accurate but higher memory footprint
Returns:

lookup_table (ndarray): Coupling lookup, indices are [rho, n1, n2]. rho_array (ndarray): Values for the radial distance considered for the lookup (starting from negative

numbers to allow for simpler cubic interpolation without distinction of cases at rho=0)
Return type:

(tuple) tuple containing
smuthi.linearsystem.particlecoupling.prepare_lookup.radial_direct_pwe_mediated_coupling_lookup_table(vacuum_wavelength, rho_max, l_max, k_is, k_parallel='default', resolution=None)

Prepare Sommerfeld integral lookup table to allow for a fast calculation of the coupling matrix by interpolation. This function is called when all particles are on the same z-position.

This function uses PVWF representation of the coupling operator and is only for the direct particle-particle lookup table for particles in close vicinity.

Parameters:
  • vacuum_wavelength (float) – Vacuum wavelength in length units
  • rho_max (float) – Maximal horizontal particle displacement for which direct coupling is evaluated via PVWFs
  • l_max (int) – Maximal multipole degree
  • k_is (complex) – Wavenumber of the scattering layer
  • k_parallel (numpy.ndarray or str) – In-plane wavenumber for Sommerfeld integrals. If ‘default’, smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float) – Spatial resolution of lookup table in length units. (default: vacuum_wavelength / 100) Smaller means more accurate but higher memory footprint
Returns:

lookup_table (ndarray): Coupling lookup, indices are [rho, n1, n2]. rho_array (ndarray): Values for the radial distance considered for the lookup (starting from negative

numbers to allow for simpler cubic interpolation without distinction of cases at rho=0)
Return type:

(tuple) tuple containing
smuthi.linearsystem.particlecoupling.prepare_lookup.size_format(b)
smuthi.linearsystem.particlecoupling.prepare_lookup.volumetric_coupling_lookup_table(vacuum_wavelength, particle_list, layer_system, k_parallel='default', resolution=None)

Prepare Sommerfeld integral lookup table to allow for a fast calculation of the coupling matrix by interpolation. This function is called when not all particles are on the same z-position.

Parameters:
  • vacuum_wavelength (float) – Vacuum wavelength in length units
  • particle_list (list) – List of particle objects
  • layer_system (smuthi.layers.LayerSystem) – Stratified medium
  • k_parallel (numpy.ndarray or str) – In-plane wavenumber for Sommerfeld integrals. If ‘default’, smuthi.fields.default_Sommerfeld_k_parallel_array
  • resolution (float) – Spatial resolution of lookup table in length units. (default: vacuum_wavelength / 100) Smaller means more accurate but higher memory footprint
Returns:

tuple containing:

w_pl (ndarray): Coupling lookup for z1 + z2, indices are [rho, z, n1, n2]. Includes layer mediated coupling. w_mn (ndarray): Coupling lookup for z1 + z2, indices are [rho, z, n1, n2]. Includes layer mediated and

direct coupling.

rho_array (ndarray): Values for the radial distance considered for the lookup (starting from negative

numbers to allow for simpler cubic interpolation without distinction of cases for lookup edges

sz_array (ndarray): Values for the sum of z-coordinates (z1 + z2) considered for the lookup dz_array (ndarray): Values for the difference of z-coordinates (z1 - z2) considered for the lookup
Return type:

(tuple)

particlecoupling.prepare_lookup_cuda

This module contains CUDA source code for the preparation of coupling matrix lookups.